function SimulatedMoments  = v10_Moments(ParametersIn);

% Parameters to be estimated

epsilon        =    ParametersIn(1);

% Parameters carried over from previous estimation




gamma = 43.14;
alpha = 6.76;
c1 = 169.16;
c2 = 41690;
b1 = 76.41;
b2 = 49.90;
z = 3.6;

gamma1 = gamma + (epsilon * 3.39);
gamma12 = gamma - (epsilon * 3.79);


%Set parameters
p_a = 1;
p_b = 0;
p_c = 0.5;

% True reporting rates
q_a = 0.900;
q_b = 0.409;
q_c = 0.507;
q_d = 0.234;

%Perceived reporting rates (March 24)
q_hat_a = 0.88;
q_hat_b = 0.59;
q_hat_c = 0.77;
q_hat_d = 0.23;

%Truncated true reporting costs
k_a =  9.690 ;
k_b =  4.065 ;
k_c = 5.196 ;
k_d = 2.054 ;

%Truncated perceived reporting costs
k_hat_a = 9.461 ;
k_hat_b = 6.149 ;
k_hat_c =  8.205 ;
k_hat_d =  2.007  ;



m = 20 ;


% Case when experience is high

% PS effort  if s=0

num2_a =  alpha*gamma12*  ((q_hat_a*((m*p_a)-k_hat_a))+b1)  *  ((q_hat_a*m*(1-p_a))+b2);
num2_b =  alpha*gamma12*  ((q_hat_b*((m*p_b)-k_hat_b))+b1)  *  ((q_hat_b*m*(1-p_b))+b2);
num2_c =  alpha*gamma12*  ((q_hat_c*((m*p_c)-k_hat_c))+b1)  *  ((q_hat_c*m*(1-p_c))+b2);
num2_d =  alpha*gamma12*  ((q_hat_d*(-k_hat_d))+b1)  *  b2;


den_a = (2*c1*c2) - ((gamma12^2)* ((q_hat_a*((m*p_a)-k_hat_a))+b1)  *  ((q_hat_a*m*(1-p_a))+b2));
den_b = (2*c1*c2) - ((gamma12^2)* ((q_hat_b*((m*p_b)-k_hat_b))+b1)  *  ((q_hat_b*m*(1-p_b))+b2));
den_c = (2*c1*c2) - ((gamma12^2)* ((q_hat_c*((m*p_c)-k_hat_c))+b1)  *  ((q_hat_c*m*(1-p_c))+b2));
den_d = (2*c1*c2) - ((gamma12^2)* ((q_hat_d*(-k_hat_d))+b1)  * b2);


e2n_a = num2_a / den_a;
e2n_b = num2_b / den_b;
e2n_c = num2_c / den_c;
e2n_d = num2_d / den_d;


% CHW effort if s=0


num1_a = alpha*((q_a*((m*p_a)-k_a))+b1)*c2;
num1_b = alpha*((q_b*((m*p_b)-k_b))+b1)*c2;
num1_c = alpha*((q_c*((m*p_c)-k_c))+b1)*c2;
num1_d = alpha*((q_d*(-k_d))+b1)*c2;


e1n_a = num1_a / den_a;
e1n_b = num1_b / den_b;
e1n_c = num1_c / den_c;
e1n_d = num1_d / den_d;


% PS effort if s>0

eta_a = (b1*z) +b2 + q_hat_a*((m* (1+ (p_a*(z-1)))-(k_hat_a*z) ));
eta_b = (b1*z) +b2 + q_hat_b*((m* (1+ (p_b*(z-1)))-(k_hat_b*z) ));
eta_c = (b1*z) +b2 + q_hat_c*((m* (1+ (p_c*(z-1)))-(k_hat_c*z) ));
eta_d = (b1*z) +b2 + q_hat_d*( -(k_hat_d*z) );


dens_a = (8*z*c1*c2)-((gamma12^2)*(eta_a^2));
dens_b = (8*z*c1*c2)-((gamma12^2)*(eta_b^2));
dens_c = (8*z*c1*c2)-((gamma12^2)*(eta_c^2));
dens_d = (8*z*c1*c2)-((gamma12^2)*(eta_d^2));


nums2_a = gamma12*alpha*(eta_a^2);
nums2_b = gamma12*alpha*(eta_b^2);
nums2_c = gamma12*alpha*(eta_c^2);
nums2_d = gamma12*alpha*(eta_d^2);

e2s_a = nums2_a / dens_a;
e2s_b = nums2_b / dens_b;
e2s_c = nums2_c / dens_c;
e2s_d = nums2_d / dens_d;


% CHW effort if s>0

e1term1_a = (4*alpha*z*c2) / dens_a;
e1term1_b = (4*alpha*z*c2) / dens_b;
e1term1_c = (4*alpha*z*c2) / dens_c;
e1term1_d = (4*alpha*z*c2) / dens_d;


e1term2_a = (q_a*(eta_a + (2*q_hat_a*z*(k_hat_a-k_a))))/ (2*q_hat_a*z);
e1term2_b = (q_b*(eta_b + (2*q_hat_b*z*(k_hat_b-k_b))))/ (2*q_hat_b*z);
e1term2_c = (q_c*(eta_c + (2*q_hat_c*z*(k_hat_c-k_c))))/ (2*q_hat_c*z);
e1term2_d = (q_d*(eta_d + (2*q_hat_d*z*(k_hat_d-k_d))))/ (2*q_hat_d*z);

e1term3_a = ( b1 * (q_hat_a-q_a)) / q_hat_a ;
e1term3_b = ( b1 * (q_hat_b-q_b)) / q_hat_b ;
e1term3_c = ( b1 * (q_hat_c-q_c)) / q_hat_c ;
e1term3_d = ( b1 * (q_hat_d-q_d)) / q_hat_d ;

e1s_a = e1term1_a * (e1term2_a + e1term3_a);
e1s_b = e1term1_b * (e1term2_b + e1term3_b);
e1s_c = e1term1_c * (e1term2_c + e1term3_c);
e1s_d = e1term1_d * (e1term2_d + e1term3_d);





% Calculate condition for s=1 
condition_a = (b2+(q_hat_a*m)-(z*b1)+(q_hat_a*k_hat_a*z)) / (q_hat_a*m*(z+1));
condition_b = (b2+(q_hat_b*m)-(z*b1)+(q_hat_b*k_hat_b*z)) / (q_hat_b*m*(z+1));
condition_c = (b2+(q_hat_c*m)-(z*b1)+(q_hat_c*k_hat_c*z)) / (q_hat_c*m*(z+1));


% Calculate output
e1_d = e1n_d;
e2_d = e2n_d;
x_d12 = alpha* e1_d + (gamma12 * e1_d * e2_d);


if p_a <= condition_a
    e1_a = e1s_a;
    e2_a = e2s_a;
    x_a12 = alpha* e1_a + (gamma12 * e1_a * e2_a);
else
    e1_a = e1n_a;
    e2_a = e2n_a;
    x_a12 = alpha* e1_a + (gamma12 * e1_a * e2_a);

end


if p_b <= condition_b
    e1_b = e1s_b;
    e2_b = e2s_b;
    x_b12 = alpha* e1_b + (gamma12 * e1_b * e2_b);

else
    e1_b = e1n_b;
    e2_b = e2n_b;
    x_b12 = alpha* e1_b + (gamma12 * e1_b * e2_b);

end

if p_c <= condition_c
    e1_c = e1s_c;
    e2_c = e2s_c;
    x_c12 = alpha* e1_c + (gamma12 * e1_c * e2_c);

else
    e1_c = e1n_c;
    e2_c = e2n_c;
    x_c12 = alpha* e1_c + (gamma12 * e1_c * e2_c);

end


% Case when experience is low



% PS effort  if s=0

num2_a =  alpha*gamma1*  ((q_hat_a*((m*p_a)-k_hat_a))+b1)  *  ((q_hat_a*m*(1-p_a))+b2);
num2_b =  alpha*gamma1*  ((q_hat_b*((m*p_b)-k_hat_b))+b1)  *  ((q_hat_b*m*(1-p_b))+b2);
num2_c =  alpha*gamma1*  ((q_hat_c*((m*p_c)-k_hat_c))+b1)  *  ((q_hat_c*m*(1-p_c))+b2);
num2_d =  alpha*gamma1*  ((q_hat_d*(-k_hat_d))+b1)  *  b2;


den_a = (2*c1*c2) - ((gamma1^2)* ((q_hat_a*((m*p_a)-k_hat_a))+b1)  *  ((q_hat_a*m*(1-p_a))+b2));
den_b = (2*c1*c2) - ((gamma1^2)* ((q_hat_b*((m*p_b)-k_hat_b))+b1)  *  ((q_hat_b*m*(1-p_b))+b2));
den_c = (2*c1*c2) - ((gamma1^2)* ((q_hat_c*((m*p_c)-k_hat_c))+b1)  *  ((q_hat_c*m*(1-p_c))+b2));
den_d = (2*c1*c2) - ((gamma1^2)* ((q_hat_d*(-k_hat_d))+b1)  * b2);


e2n_a = num2_a / den_a;
e2n_b = num2_b / den_b;
e2n_c = num2_c / den_c;
e2n_d = num2_d / den_d;


% CHW effort if s=0


num1_a = alpha*((q_a*((m*p_a)-k_a))+b1)*c2;
num1_b = alpha*((q_b*((m*p_b)-k_b))+b1)*c2;
num1_c = alpha*((q_c*((m*p_c)-k_c))+b1)*c2;
num1_d = alpha*((q_d*(-k_d))+b1)*c2;


e1n_a = num1_a / den_a;
e1n_b = num1_b / den_b;
e1n_c = num1_c / den_c;
e1n_d = num1_d / den_d;


% PS effort if s>0

eta_a = (b1*z) +b2 + q_hat_a*((m* (1+ (p_a*(z-1)))-(k_hat_a*z) ));
eta_b = (b1*z) +b2 + q_hat_b*((m* (1+ (p_b*(z-1)))-(k_hat_b*z) ));
eta_c = (b1*z) +b2 + q_hat_c*((m* (1+ (p_c*(z-1)))-(k_hat_c*z) ));
eta_d = (b1*z) +b2 + q_hat_d*( -(k_hat_d*z) );


dens_a = (8*z*c1*c2)-((gamma1^2)*(eta_a^2));
dens_b = (8*z*c1*c2)-((gamma1^2)*(eta_b^2));
dens_c = (8*z*c1*c2)-((gamma1^2)*(eta_c^2));
dens_d = (8*z*c1*c2)-((gamma1^2)*(eta_d^2));


nums2_a = gamma1*alpha*(eta_a^2);
nums2_b = gamma1*alpha*(eta_b^2);
nums2_c = gamma1*alpha*(eta_c^2);
nums2_d = gamma1*alpha*(eta_d^2);

e2s_a = nums2_a / dens_a;
e2s_b = nums2_b / dens_b;
e2s_c = nums2_c / dens_c;
e2s_d = nums2_d / dens_d;


% CHW effort if s>0

e1term1_a = (4*alpha*z*c2) / dens_a;
e1term1_b = (4*alpha*z*c2) / dens_b;
e1term1_c = (4*alpha*z*c2) / dens_c;
e1term1_d = (4*alpha*z*c2) / dens_d;


e1term2_a = (q_a*(eta_a + (2*q_hat_a*z*(k_hat_a-k_a))))/ (2*q_hat_a*z);
e1term2_b = (q_b*(eta_b + (2*q_hat_b*z*(k_hat_b-k_b))))/ (2*q_hat_b*z);
e1term2_c = (q_c*(eta_c + (2*q_hat_c*z*(k_hat_c-k_c))))/ (2*q_hat_c*z);
e1term2_d = (q_d*(eta_d + (2*q_hat_d*z*(k_hat_d-k_d))))/ (2*q_hat_d*z);

e1term3_a = ( b1 * (q_hat_a-q_a)) / q_hat_a ;
e1term3_b = ( b1 * (q_hat_b-q_b)) / q_hat_b ;
e1term3_c = ( b1 * (q_hat_c-q_c)) / q_hat_c ;
e1term3_d = ( b1 * (q_hat_d-q_d)) / q_hat_d ;

e1s_a = e1term1_a * (e1term2_a + e1term3_a);
e1s_b = e1term1_b * (e1term2_b + e1term3_b);
e1s_c = e1term1_c * (e1term2_c + e1term3_c);
e1s_d = e1term1_d * (e1term2_d + e1term3_d);




% Calculate condition for s=1 
condition_a = (b2+(q_hat_a*m)-(z*b1)+(q_hat_a*k_hat_a*z)) / (q_hat_a*m*(z+1));
condition_b = (b2+(q_hat_b*m)-(z*b1)+(q_hat_b*k_hat_b*z)) / (q_hat_b*m*(z+1));
condition_c = (b2+(q_hat_c*m)-(z*b1)+(q_hat_c*k_hat_c*z)) / (q_hat_c*m*(z+1));


% Calculate output
e1_d = e1n_d;
e2_d = e2n_d;
x_d1 = alpha* e1_d + (gamma1 * e1_d * e2_d);


if p_a <= condition_a
    e1_a = e1s_a;
    e2_a = e2s_a;
    x_a1 = alpha* e1_a + (gamma1 * e1_a * e2_a);
else
    e1_a = e1n_a;
    e2_a = e2n_a;
    x_a1 = alpha* e1_a + (gamma1 * e1_a * e2_a);

end


if p_b <= condition_b
    e1_b = e1s_b;
    e2_b = e2s_b;
    x_b1 = alpha* e1_b + (gamma1 * e1_b * e2_b);

else
    e1_b = e1n_b;
    e2_b = e2n_b;
    x_b1 = alpha* e1_b + (gamma1 * e1_b * e2_b);

end

if p_c <= condition_c
    e1_c = e1s_c;
    e2_c = e2s_c;
    x_c1 = alpha* e1_c + (gamma1 * e1_c * e2_c);

else
    e1_c = e1n_c;
    e2_c = e2n_c;
    x_c1 = alpha* e1_c + (gamma1 * e1_c * e2_c);

end


% Calculate moments

diff_a = ((x_a1 - x_d1) - (x_a12 - x_d12))/7.18;
diff_b = ((x_b1 - x_d1) - (x_b12 - x_d12))/7.18;
diff_c = ((x_c1 - x_d1) - (x_c12 - x_d12))/7.18;

SimulatedMoments    = [diff_a; diff_b; diff_c];

return
